Does interfertility decrease with phylogenetic distance in Schiedea?
Binary pollen fertility data
Data from Weller et al. 2001. F1 hybrids that made greater or less than 50 percent fertile pollen.
fert <- read.table(file="data/fertility.csv", sep="\t", header=T)
fert$sp1 <- toupper(fert$sp1)
fert$sp2 <- toupper(fert$sp2)
fert$g50 <- fert$g50+1
fspecies <- data.frame(name=levels(factor(c(fert$sp1,fert$sp2))))
g <- graph_from_data_frame(fert, directed=F, vertices=fspecies)
g_edges = get.edgelist(g)Clustering
Phylogeny on the right, clustering by network of fertile pollen on left
#clustering by network
adj <- get.adjacency(g, attr="g50")
adj[adj==2] <- 0.5
g.d <- graph_from_adjacency_matrix(adj, weighted=T, mode="undirected")
ceb <- cluster_edge_betweenness(g.d)
#trim phlyo
sch.f <- sch
sch.f$tip.label <- toupper(paste("",substr(sch.f$tip.label,1,4),sep=""))
sch.f<-drop.tip(sch.f,sch.f$tip.label[-match(fspecies$name, sch$tip.label)])
#cophenetic plot
library(phytools)
ceb_tree <- cophylo(as_phylo(ceb), sch.f)Rotating nodes to optimize matching...
Done.plot(ceb_tree)
Arcplot
F1 hybrids that made greater (green) or less than (red) 50 percent fertile pollen.
#phylo-arcplot
par(mfrow=c(1,2))
par(mar=c(1.2,0,1.2,0))
plot(sch.f)
par(mar=c(0,0,0,0))
arcplot(g_edges, show.labels=F, horizontal=F, ordering=sch.f$tip.label, col.arcs=fert$g50+1)
Fertility vs. phylogenetic distance
adj <- as.matrix(get.adjacency(g, attr="g50"))
adj[upper.tri(adj, diag=T)]<-NA
adj[adj==0] <- NA
phydist <- cophenetic.phylo(sch.f)
phydist[upper.tri(phydist, diag=T)]<-NA
pf <- na.omit(data.frame(pd=c(phydist), ft=c(adj-1)))
cor(pf$ft, pf$pd, use="complete.obs")[1] -0.3740546cor.test(pf$ft, pf$pd, use="complete.obs")
Pearson's product-moment correlation
data: pf$ft and pf$pd
t = -2.9085, df = 52, p-value = 0.005331
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.5833851 -0.1181252
sample estimates:
cor
-0.3740546 Logistic regression
#Logistic regression
model <- glm(ft~pd, family=binomial(link='logit'), data=pf)
summary(model)
Call:
glm(formula = ft ~ pd, family = binomial(link = "logit"), data = pf)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.546 1.689 2.691 0.00712 **
pd -177.869 70.976 -2.506 0.01221 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 71.188 on 53 degrees of freedom
Residual deviance: 62.457 on 52 degrees of freedom
AIC: 66.457
Number of Fisher Scoring iterations: 4confint(model) 2.5 % 97.5 %
(Intercept) 1.660361 8.395113
pd -337.364722 -54.308760wald.test(b = coef(model), Sigma = vcov(model), Terms = 2)Wald test:
----------
Chi-squared test:
X2 = 6.3, df = 1, P(> X2) = 0.012Logistic regression plot
plot <- ggplot(pf, aes(pd, ft)) + stat_smooth(method="glm", method.args=list(family="binomial"), se=T, color="black") + scale_y_continuous(breaks=c(0, 1), labels=c("<50%", ">50%")) + scale_x_continuous() + xlab("Phylogenetic distance") + ylab("Hybrid pollen fertility") + scale_fill_gradientn(colors=c("white","black")) + theme_bw()
plot + geom_bin2d()
Continuous pollen fertility data
crosses <- read.table("data/hpol_crosses.csv", sep="\t", header=T, colClasses = list(cross="factor",crosstype="factor"))
str(crosses)'data.frame': 120 obs. of 11 variables:
$ cross : Factor w/ 120 levels "108","11","116",..: 13 24 105 107 118 2 4 9 20 23 ...
$ mompop : int 847 847 847 847 847 847 847 847 847 847 ...
$ momid : chr "B" "B" "B" "B" ...
$ dadpop : int 794 794 794 794 794 842 842 842 847 791 ...
$ dadid : chr "3" "3" "3" "3" ...
$ momsp : chr "adam" "adam" "adam" "adam" ...
$ dadsp : chr "hook" "hook" "hook" "hook" ...
$ mombs : chr "d" "d" "d" "d" ...
$ dadbs : chr "h" "h" "h" "h" ...
$ crosstype: Factor w/ 4 levels "1","2","3","4": 3 3 3 3 3 4 4 4 4 4 ...
$ notes : chr "" "" "" "" ...pol <- read.table("data/hpol.csv", sep="\t", header=T, colClasses = list(cross="factor", id="factor",crosstype="factor"))
pol$vp <- pol$viable/200
pol$mombs <- factor(ifelse(pol$crosstype %in% c("1","2"), "H", ifelse(pol$crosstype=="z", "z", "D")))
pol$dadbs <- factor(ifelse(pol$crosstype %in% c("1","3"), "H", ifelse(pol$crosstype=="z", "z", "D")))
pol$momgrp <- with(pol, factor(paste(momclade,momisle,mombs)))
pol$dadgrp <- with(pol, factor(paste(dadclade,dadisle,dadbs)))
pol$momsp <- crosses$momsp[match(pol$cross, crosses$cross)]
pol$dadsp <- crosses$dadsp[match(pol$cross, crosses$cross)]
str(pol)'data.frame': 2239 obs. of 16 variables:
$ cross : Factor w/ 113 levels "108","11","116",..: 99 99 99 99 99 10 10 10 10 10 ...
$ id : Factor w/ 38 levels "1","10","11",..: 25 25 25 25 25 11 11 11 11 11 ...
$ viable : int 122 165 151 166 173 189 191 183 184 179 ...
$ crosstype: Factor w/ 5 levels "-9","1","2","3",..: 5 5 5 5 5 5 5 5 5 5 ...
$ date : int 1 1 1 1 1 -9 -9 -9 -9 -9 ...
$ momisle : chr "O" "O" "O" "O" ...
$ dadisle : chr "O" "O" "O" "O" ...
$ momclade : chr "2" "2" "2" "2" ...
$ dadclade : chr "4" "4" "4" "4" ...
$ vp : num 0.61 0.825 0.755 0.83 0.865 0.945 0.955 0.915 0.92 0.895 ...
$ mombs : Factor w/ 2 levels "D","H": 1 1 1 1 1 1 1 1 1 1 ...
$ dadbs : Factor w/ 2 levels "D","H": 1 1 1 1 1 1 1 1 1 1 ...
$ momgrp : Factor w/ 13 levels "1 K H","1 N H",..: 5 5 5 5 5 11 11 11 11 11 ...
$ dadgrp : Factor w/ 16 levels "1 K H","1 N H",..: 13 13 13 13 13 5 5 5 5 5 ...
$ momsp : chr "adam" "adam" "adam" "adam" ...
$ dadsp : chr "keal" "keal" "keal" "keal" ...heatmap(with(pol, table(momsp,dadsp)), col=inferno(512), scale="none", Colv=NA, Rowv=NA)
#polmeans <- aggregate(vp~cross*momgrp*momclade*momisle*mombs*dadgrp*dadclade*dadisle*dadbs, data=pol, FUN=function(x) c(mean(x), sd(x), length(x)))
polmeans <- aggregate(vp~momsp*dadsp, data=pol, FUN=function(x) c(median(x), quantile(x,0.25), quantile(x,0.75), length(x))) #mean is median!!!
polmeans <- do.call(data.frame, polmeans)
colnames(polmeans)[3:6] <- c("mean", "q25", "q75", "N")
#write.csv(polmeans, "results/hpol_means.csv", row.names=F)
#alllvl <- sort(unique(c(levels(polmeans$momgrp), levels(polmeans$dadgrp))))
#grps <- cbind(alllvl, read.table(text=alllvl, sep=" ", col.names=c("clade","isle","bs")))
species <- unique(c(pol$momsp,pol$dadsp))
grps <- data.frame(sp=species)
library(reshape2)
poltbl <- dcast(polmeans, momsp ~ dadsp, value.var="mean", fun.aggregate = mean)
rownames(poltbl) <- poltbl$momsp
poltbl$momsp <- NULL
poltbl2 <- as.matrix(poltbl)
poltbl2[is.nan(poltbl2)] <- 0ggplot(pol, aes(x=momsp, color=dadsp, y=vp)) + geom_boxplot()
ggplot(pol, aes(x=momclade, color=dadclade, y=vp)) + geom_boxplot() + geom_point()
ggplot(pol, aes(x=momisle, color=dadisle, y=vp)) + geom_boxplot() + geom_point()
ggplot(pol, aes(x=momclade, color=dadclade, y=vp)) + geom_boxplot() + geom_point()
ggplot(pol, aes(x=reorder(cross, viable, FUN=mean), y=vp, color=crosstype)) + geom_boxplot() + geom_point()
ggplot(pol, aes(x=crosstype, y=vp, color=dadclade, fill=momclade)) + geom_boxplot()
#heatmap(poltbl2, col=viridis(512), scale="none", Rowv = NA, Colv=NA)
#ggplot(polmeans, aes(x=factor(momsp, levels=alllvl), y=factor(dadgrp, levels=alllvl), fill=mean)) + geom_tile() + scale_x_discrete(drop=F) + scale_y_discrete(drop=F) + scale_fill_viridis_c()
#ggplot(polmeans, aes(x=momclade, y=dadclade, fill=mean)) + geom_tile() + scale_fill_viridis_c()
ggplot(polmeans, aes(x=momsp, y=dadsp, fill=mean)) + geom_tile() + scale_fill_viridis_c()
hist(pol$viable, breaks=40)
colfunc <- viridisLite::plasma
g <- graph_from_data_frame(polmeans[,c("momsp","dadsp","mean","q25","q75","N")], directed=T, vertices=grps)
lgw <- layout_in_circle(g)
#plot.igraph(g, layout=layout_in_circle(g, order=order(grps$clade)), edge.color=colfunc(100)[round(E(g)$mean*100)], vertex.color=rainbow(10)[grps$clade], vertex.label.color=rainbow(10)[as.integer(grps$isle)+6], vertex.label.cex=2)
plot.igraph(g, layout=layout_in_circle(g))
gu <- as.undirected(g, mode="collapse", edge.attr.comb = "mean")
adj <- get.adjacency(gu, attr="mean")
g.w <- graph_from_adjacency_matrix(adj, weighted=T, mode="undirected")
#plot(g.w, layout_with_fr(g.w), edge.color=colfunc(100)[round(E(g)$mean*100)], edge.width=2, vertex.color=rainbow(10)[grps$clade], vertex.label.color=rainbow(10)[as.integer(grps$isle)+6], vertex.label.cex=2)
plot(g.w, layout=layout_with_fr(g.w), edge.color=colfunc(100)[round(E(g)$mean*100)], edge.width=2)
#dendPlot(cle, mode="phylo", use.edge.length = T, main="leading eigen")
plot(as.dendrogram(cluster_leading_eigen(gu)))
g_edges = get.edgelist(g)
arcplot(g_edges, show.labels=T, horizontal=F, col.arcs=colfunc(100)[round(polmeans$mean*100)], sorted=T)
Clustering
ceb <- cluster_leading_eigen(gu)
#trim phlyo
sch.f <- sch
sch.f$tip.label <- tolower(paste("",substr(sch.f$tip.label,1,4),sep=""))
sch.f<-drop.tip(sch.f,sch.f$tip.label[-match(species, sch.f$tip.label)])
species1f <- species1
rownames(species1f) <- tolower(rownames(species1f))
species1f <- species1f[sch.f$tip.label,]
#BScols <- setNames(c("black", "navy", "lightblue", "white", "white", "white"),1:6)
BScols <- setNames(RColorBrewer::brewer.pal(6,"Set2")[c(6,4,2,5,1,3)],1:6)
#intraspecies <- setNames(polmeans$mean[match(paste(sch.f$tip.label,sch.f$tip.label), paste(polmeans$momsp,polmeans$dadsp))], sch.f$tip.label)
intraspecies <- setNames(polmeans$mean[match(paste(names(V(g)),names(V(g))), paste(polmeans$momsp,polmeans$dadsp))], names(V(g)))
#cophenetic plot
ceb_tree <- cophylo(as_phylo(ceb), sch.f)Rotating nodes to optimize matching...
Done.plot(ceb_tree)
#phylo-arcplot
#png("plots/arcplot.png", height=1500, width=3000, pointsize=56)
par(mfrow=c(1,3) , bg="grey30")
par(mar=c(1.5,0,1.5,0))
plot(sch.f, align.tip.label=T, tip.color=BScols[as.integer(species1f$Breeding.System)], edge.width=8, edge.color="grey80", cex=1.5, lwd=8)
add.scale.bar(y=0.002, x=0.002, lcol="white",col="grey80")
par(mar=c(0,0,0,0))
arcplot(g_edges, show.labels=F, horizontal=F, ordering=sch.f$tip.label, col.arcs=plasma(512, direction=-1)[round(polmeans$mean*512)], sorted=T, show.nodes=T, lwd.arcs=6, col.nodes = plasma(512, direction=-1)[round(intraspecies*512)], pch.nodes = 19, cex.nodes=1)
plot.new()
#legend(x=-0.05,y=1, legend = rev(levels(species1f$Breeding.System)), col = rev(BScols), pch = 19, cex=1.2, pt.cex=2.5, bty="n", xpd=T, text.col= rev(BScols))
library(plotrix)
color.legend(0,0.55,1,0.65,legend=seq(0,1,0.2),rect.col=plasma(512, direction=-1),cex=1, col="grey80", align="rt")
text(0.2, 0.7, "Pollen viability", col="grey80", cex=1.5)
#dev.off()Fertility vs. phylogenetic distance
phydist <- cophenetic.phylo(sch.f)
#phydist[upper.tri(phydist, diag=T)]<-NA
adj <- as.matrix(as_adj(g, attr="mean"))
#diag(adj) <- NA
#adj[upper.tri(adj, diag=T)]<-NA
adj[adj==0] <- NA
adj <- adj[sch.f$tip.label,sch.f$tip.label]
adj.N <- as.matrix(as_adj(g, attr="N"))
adj.N <- adj.N[sch.f$tip.label,sch.f$tip.label]
adj.q25 <- as.matrix(as_adj(g, attr="q25"))
adj.q25 <- adj.q25[sch.f$tip.label,sch.f$tip.label]
adj.q75 <- as.matrix(as_adj(g, attr="q75"))
adj.q75 <- adj.q75[sch.f$tip.label,sch.f$tip.label]
comb <-function(x,y) {mapply(function(q,r) paste(sort(c(q,r)), collapse="-"), x, y)}
labels <- outer(sch.f$tip.label, sch.f$tip.label, comb)
pf <- na.omit(data.frame(pd=c(phydist), ft=c(adj), N=c(adj.N), q25=c(adj.q25), q75=c(adj.q75),crossname= c(outer(sch.f$tip.label, sch.f$tip.label, paste, sep=" x ")), crossname.sort=c(labels)))
pf$q75[pf$q75==1]<- NA # safe for logit transform
cor(pf$ft, pf$pd, use="complete.obs")[1] -0.7844253cor.test(pf$ft, pf$pd, use="complete.obs")
Pearson's product-moment correlation
data: pf$ft and pf$pd
t = -11.026, df = 76, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.8573058 -0.6807202
sample estimates:
cor
-0.7844253 pdft <- ggplot(pf, aes(x=pd, y=ft, weight=200*N)) +
geom_line(aes(group=crossname.sort), col="grey70", linewidth=2, alpha=0.5) +
theme_grey() +
xlab("Phylogenetic distance (mutations/site)") +
ylab("Pollen viability")
pdft + geom_point() + geom_linerange(aes(ymin=q25, ymax=q75)) + ylim(c(0,1))
pdft + scale_y_continuous(trans="logit", breaks=c(0.01,seq(0,1,by=0.05), 0.995)) + geom_text(aes(label=crossname), size=3) # trans_breaks(function(p) { log(p/(1-p)) }, function(x) { 1/(1+exp(-x)) }, n=10)
pdft + geom_smooth(method="glm", method.args=list(family="binomial"), color="grey70", fill="grey80") + geom_text(aes(label=crossname, color=N)) + coord_trans(y="logit") + scale_y_continuous(breaks=c(0.01,seq(0,1,by=0.05), 0.994)) + scale_color_gradient(trans="sqrt", low="darkblue", high="green2")
#scale_y_continuous(breaks=seq(0,1, by=0.05), limits=c(0,1), expand=c(0,0))
#ggsave("pdft.png", height=9, width=9)Arcplot
#phylo-arcplot
#TODO: arrange arcs front-to-back by phydist?
library(gridBase)
library(grid)
#png("plots/arcplot.png", height=1500, width=4000, pointsize=20, res=300)
par(mfrow=c(1,4) , bg="grey30")
par(mar=c(1,0,1,0))
plot(sch.f, align.tip.label=T, tip.color=BScols[as.integer(species1f$Breeding.System)], edge.width=3, edge.color="grey80", cex=1, lwd=4, label.offset=0.0001)
add.scale.bar(y=0.6, x=0, lcol="grey80",col="grey80", cex=1, lwd=3)
par(mar=c(0,0,0,0))
arcplot(g_edges, show.labels=F, horizontal=F, ordering=sch.f$tip.label, col.arcs=plasma(512, direction=-1)[round(polmeans$mean*512)], sorted=T, show.nodes=T, lwd.arcs=3, col.nodes = plasma(512, direction=-1)[round(intraspecies[order(names(intraspecies))]*512)], pch.nodes = 19, cex.nodes=1)
plot.new()
#legend(x=-0.05,y=1, legend = rev(levels(species1f$Breeding.System)), col = rev(BScols), pch = NA, cex=1, pt.cex=2.5, bty="n", xpd=T, text.col= rev(BScols))
library(plotrix)
color.legend(0.1,0.05,0.2,0.55,legend=seq(0,1,0.2),rect.col=plasma(511, direction=-1),cex=0.45, col="grey80", align="rt", gradient="y")
text(0.26, 0.6, "Pollen viability", col="grey80", cex=1)
plot.new()
vps <- baseViewports()
pushViewport(vps$figure)
vp1 <-plotViewport(c(0,0,0,0))
print(pdft +
geom_smooth(method="glm", method.args=list(family="binomial"), color="grey70", fill="grey80") +
geom_linerange(aes(ymin=q25, ymax=q75)) +
geom_point(aes(color=ft), size=2) +
coord_trans(y="logit") +
scale_y_continuous(breaks=c(0.01,seq(0,1,by=0.1), 0.99)) +
scale_color_viridis_c(option="plasma", direction=-1) +
theme_classic() +
theme(plot.background = element_rect(color = "grey30", fill = "grey30"), panel.background = element_rect(fill = "grey30", color = NA), legend.background = element_rect(color = NA, fill = "grey30"), axis.text = element_text(color = "grey80"), axis.title = element_text(color = "grey80")) +
guides(color=F), vp = vp1)
#dev.off()Binomial regression on means
# model <- glm(ft~pd, family="quasibinomial", weights=N*200, data=pf)
# summary(model)
#plot(model)
library(glmmTMB)
betabin.tmb <- glmmTMB(round(cbind(ft*N*200,(1-ft)*N*200))~pd, data=pf, family="betabinomial")
summary(betabin.tmb) Family: betabinomial ( logit )
Formula: round(cbind(ft * N * 200, (1 - ft) * N * 200)) ~ pd
Data: pf
AIC BIC logLik deviance df.resid
1220.9 1227.9 -607.4 1214.9 75
Dispersion parameter for betabinomial family (): 3.63
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.9620 0.3353 8.833 < 2e-16 ***
pd -134.2268 16.5402 -8.115 4.85e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Binomial regression on raw data
phydist.m <- melt(phydist)
polmeans$pd <- phydist.m$value[match(paste(polmeans$momsp,polmeans$dadsp), paste(phydist.m$Var1,phydist.m$Var2))]
pol$pd <- phydist.m$value[match(paste(pol$momsp,pol$dadsp), paste(phydist.m$Var1,phydist.m$Var2))]
pol$tot <- 200
pol$inviable <- pol$tot-pol$viable
ggplot(pol, aes(y=vp,x=pd)) + geom_point()
# model <- glm(cbind(viable, inviable)~pd, family="quasibinomial", data=pol)
# summary(model)
# #plot(model)
# library(jtools)
# effect_plot(model, pred = pd, interval = TRUE, data=pol) + geom_point(data=pol, aes(pd,vp), color="blue")
# coef(model)
library(glmmTMB)
# plot.glmmTMB <- function(x, ...) { qplot(x$fitted, x$residuals, xlab="Fitted values", ylab="Residuals", ...) + geom_smooth()}
#
# betabin.tmb <- glmmTMB(cbind(viable,inviable)~pd, data=pol, family="betabinomial")
# summary(betabin.tmb)
#plot(betabin.tmb)
# plot(fitted(betabin.tmb)~pol$pd)
# points(polmeans$mean~polmeans$pd, col="red")
#coef(betabin.tmb)
betabin.tmb.momdadsp <- glmmTMB(cbind(viable,inviable)~pd + (1|momsp) + (1|dadsp), data=pol, family="betabinomial")
#load("betabin.tmb.momdadsp.Rdata")
summary(betabin.tmb.momdadsp) Family: betabinomial ( logit )
Formula: cbind(viable, inviable) ~ pd + (1 | momsp) + (1 | dadsp)
Data: pol
AIC BIC logLik deviance df.resid
20209.5 20238.1 -10099.8 20199.5 2234
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
momsp (Intercept) 0.8521 0.9231
dadsp (Intercept) 0.5991 0.7740
Number of obs: 2239, groups: momsp, 15; dadsp, 16
Dispersion parameter for betabinomial family (): 3.82
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.6064 0.3249 4.943 7.67e-07 ***
pd -81.4044 4.9036 -16.601 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#plot(betabin.tmb.momdadsp)
#save(betabin.tmb.momdadsp, file="betabin.tmb.momdadsp.Rdata")
######Plots############
#Looks like the random effects are overfitting :( GLM does the best
#pol$fitted <- fitted(betabin.tmb)
pol$fitted <- fitted(betabin.tmb.momdadsp)
#pol$fitted <- fitted(model)
polmeans$fitmean <- aggregate(fitted~momsp*dadsp, data=pol, FUN=mean)$fitted
ggplot(polmeans, aes(pd, fitmean)) +
geom_smooth(aes(y=mean), method="glm", method.args=list(family="binomial"), color="red", fill="pink") +
geom_smooth(method="glm", method.args=list(family="binomial"), color="black") +
geom_point() + geom_point(aes(y=mean), color="red") +
geom_segment(aes(xend=pd, yend=mean)) + geom_text(aes(y=mean, label=paste(momsp,dadsp,sep=" x ")), color="blue")
####Random effects#####
momsp.ef <- data.frame(sp=rownames(ranef(betabin.tmb.momdadsp)$cond$momsp), ef=ranef(betabin.tmb.momdadsp)$cond$momsp[[1]])#, sd=ranef(betabin.tmb.momdadsp, sd=T)$cond$momsp[[1]])
momsp.ef$sp <- reorder(momsp.ef$sp, momsp.ef$ef)
ggplot(momsp.ef, aes(sp,ef)) + geom_point() +# geom_linerange(aes(ymin=ef-sd, ymax=ef+sd)) +
ggtitle("Maternal random effects") + xlab("Species") + ylab("Effect")
dadsp.ef <- data.frame(sp=rownames(ranef(betabin.tmb.momdadsp)$cond$dadsp), ef=ranef(betabin.tmb.momdadsp)$cond$dadsp[[1]])#, sd=ranef(betabin.tmb.momdadsp, sd=T)$cond$dadsp[[1]])
dadsp.ef$sp <- reorder(dadsp.ef$sp, dadsp.ef$ef)
ggplot(dadsp.ef, aes(sp,ef)) + geom_point() +# geom_linerange(aes(ymin=ef-sd, ymax=ef+sd)) +
ggtitle("Paternal random effects") + xlab("Species") + ylab("Effect")